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The recently developped "internal" Density Functional Theory provides an existence theorem 
for a local potential that contains the center-of-mass correlations effects. The knowledge of the 
corresponding energy functional would provide a much cheaper way than projection techniques to 
treat those correlations. The aim of this article is to construct such a functional. We propose a well 
founded method, suitable for Fermions as well as for Bosons, which does not necessit to include any 
free parameter. 

I. INTRODUCTION. 

One of the most obvious symmetry of isolated self-bound systems (as atomic nuclei, Helium droplets or molecular 
systems where nuclei are treated explicitely) is translational invariance. Translational invariancc of self-bound systems 
Hamiltonians ensures Galilean invariance of the wave function, so that the center-of-mass (cm.) properties can be 
separated from the "internal" ones (that are always deduced by experimentalists). As a consequence, one laboratory 
coordinate is redundant for the description of the internal properties, which gives rise to cm. correlations. 

A numerically manageable and succesfull way to describe self-bound systems is to use mean-field like calculations 
with effective interactions. The corresponding equations are often justified starting from the Hartrce-Fock (HF) 
framework, which sacrifices by construction " Galilean invariance for the sake of the Pauli principle" , to quote Ref. [l[ . 
As a consequence, the cm. correlations (or equivalently the redundant coordinate problem) are (is) not treated 
correctly. This introduces a spurious coupling between the internal properties and the cm. motion in a HF framework, 
that affects the energy and other observables 0. 

A way to overcome this problem in the stationary case is to perform projected HF, where projection before variation 
on cm. momentum restores translational or Galilean invariance of the wave-function. Peierls and Yoccoz proposed 
a (single) projection method to restore translational invariance Q. Later, Peierls and Thouless proposed a (double) 
projection method to restore the more fundamental Galilean invarianve Q. To our knowledge, every numerical 
calculations that treat the cm. correlations by projection before variation have been done within the Peirls and 
Yoccoz method (B-Q, thus not restoring the full Galilean invariance. Moreover, the prices to pay are to abandon the 
independent-particle description and a large numerical cost [H 0; S E3]- Indeed, projection techniques require "an 
order of magnitude more computing time than the underlying mean- field like calculations" , to quote Ref. Q , which is 
discriminatory for the description of intermediary size systems. This led to the development of various approximate 
methods to treat the cm. correlations, see Ref. [TTJ] for an overview. For instance, a common method is to add a 
— < 2 ^ iN > term in the energy functional (more details will be given in mi C|) . But the successes of those methods 
are not systematic and the approximations done not perfectly justified. 

In the time-dependent case, the spurious cm. motion problem remains [l2l. [l3j. but the situation is more tricky 
because then the projected HF method becomes unmanageable even for very small self-bound systems It thus 
remains an open question to develop a rigorous and numerically cheap scheme to treat the cm. correlations, which 
would go beyond standard approximations and would remain usable in the time-dependent case. 

So far, the search for such a scheme has not been studied extensively, maybe because it is sometimes thought that 
the cm. correlations problem concerns only very small self-bound systems. But the cm correlations can have a non 
negligible effect even for intermediary size systems. For instance, it has been shown that cm. correlations are not 
negligible for all nuclei above 16 This reinforces the necessity to develop a numerically manageable method 

to treat those correlations. 

A rigorous alternative and a priori numerically much less costly way to take into account the cm. correlations has 
been opened by the recently developped "internal" Density Functional Theory (DFT) and Kohn-Sham (KS) scheme 
[lB-[I3|- Differing from standard DFT [l8i - l2lj |. it is formulated in the cm frame of a self-bound system and proves that 
the cm. correlations can be included in the energy functional and thus in a local KS potential p^l-[l7}. In addition 
to the fact that it gives a much more fundamental justification than the HF framework to the use of mean-field like 
calculations with effective interactions for the description of self-bound systems, it shows that there would be no 
need for a cm. projection if the ultimate functional was known. Internal DFT gives an existence theorem but not a 
constructive method. The aim of the present article is to propose such a method. 
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The article is organized as follows. In section II is provided a quick recall of the internal DFT formalism and 
underlined the limitations of the commonly used methods to treat the cm. correlations. In section III is developped 
a new general form for a local cm. correlations potential, whithout including any free parameter. Finally, in section 
IV are given preliminary but convincing numerical results on various model systems. 



II. INTERNAL DFT AND THE CM. CORRELATIONS FUNCTIONAL. 



A. Quick recall of the internal DFT formalism. 

We start from a self-bound system composed by N identical particles (Fermions or Bosons) and follow the consider- 
ations of Ref. [HI . The coordinates of the particles in any chosen inertial frame of reference (such as the laboratory) 
are denoted {r^}. The cm. coordinate of the system is denoted 



1 N 
N £-< J 



The system is described by the following translationally invariant iV-body Hamiltonian: 



TV N 

£ 

,3 = 1 



H = £ £ + £ u ^ + £" int (* - R ) > (!) 

i,j = l i=l 



composed by the usual kinetic energy term, a 2-body potential u which describes the Fcrmion-Fcrmion (or Boson- 
Boson) interaction (generalization of the following considerations to 3-body. . . interactions is straightforward) and an 
arbitrary translationally invariant potential v lnt . This last potential is an "internal" potential, i.e. it is defined in 
the cm. frame an acts only on the internal properties. Of course it is zero in the purely isolated self-bound case. 
Nevertheless, its form is suitable to model internal effects of fields used in experiments (polarization potentials. . . ) [jq . 

We introduce the Jacobi coordinates £ a defined as £i = r-2 — i"i, £2 = — r2 \ ri , ■ ■ ■ , £jv-i = j^zj ( r N — R)- This 
permits us to separate the Hamiltonian ([T]) into H = Hqai + Hint-, where Hcm = — (^ 2 /2A/)Ar, (with M = Nm 
being the total mass) is a 1-body Hamiltonian dcrcribing the cm. motion and acting in the R space only, and Hint 
is a (N — 1) body-Hamiltonian describing the internal properties and acting in the space only: 



N-l 9 



H in t= £^ + tf(£ 1 ,...,& v _ 1 ) + V^(£ 1 ,...,6r_ 1 ). (2) 

Zfl a 

a—l 

Hi n t contains the interaction u and the potential v mt , because they can be rewritten as functions of the {£ Q } only 
[denoted respectively f/(£i! . . . ,£n-i) and V int (£i, . . . , ^jv-i)]j an d the internal kinetic energy, which is expressed in 
terms of the conjugate momentum r a of £ Q and the reduced masses [i a = m-^-^. As [H C m, H in t] = 0, the eigenstate 
i\) of H can be build as a product of the form 

V>(ri, . . . , v N ) = T(R) Vint(6, ■ • ■ , Zn-i), (3) 
where T and ipi nt arc defined by the equations 

Hintlpint = Ei nt 1pi n t ■ (5) 

r(R) is the cm. wavefunction, that describes the motion of the isolated system as a whole in any inertial frame 
of reference. Since T is solution of the free Schrddingcr equation, T(R) should be an arbitrary stationary plane 
wave, i.e. infinitely spread and not normalizable. This leads to derealization of R and arbitrary cm. energy 
Ecm = ft K /(2M). This does not correspond to experimental situations, where the system is no longer isolated: 
interactions with other systems of the experimental apparatus localize the cm.. But the formal decoupling between 
the cm. motion and internal properties permits us to "let" the cm. motion to the choice of experimental conditions, 
internal properties described by ipint being fully comparable to the experimental ones. Note that ipi n t is by definition 
always normalizable for the ground state of a self-bound system. The internal density associated to ipint is [HI, l22l . [23j 

/ N \ 3 f 2 
Pint(r) = N [ N _ -J d£i---d£ N -2 |V>int(£l,---,6v-2, w=l)\ 
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= N fdr 1 ---dT N 5(RM int (r 1 ,...,r N )\ 2 6(r-(r i -K)). (6) 

The density Pi n t( r ) is normalized to N and r is defined in the cm. frame (see the delta relation in the previous 
equation) . Note that, even if ipi n t can be written as a function of the (N — 1) Jacobi coordinates only, it can also 
be written as a function of the N coordinates , but one of them would be redundant [HI > which is expressed by the 
<5(R) in the previous equation. 

The stationary internal DFT theorem, demonstrated by various ways in Refs. [lH HlJ [26j], states that for a non- 
degenerate ground state and for given Fermion or Boson kind of particles (i.e. a given interaction u), ipint can be 
expressed as a unique functional of pmti i- e - ^intlpint]- As a consqucncc, the ground state internal energy of a 
self-bound system Em — (i>int[pint]\Hi n t\'<Pint[pint]) can also be expressed as a unique functional of pint- 

A practical way to compute p int is given by the internal KS scheme, developped in Ref. (l5j . To set up this scheme, 
we assume that there exists, in the c.m. frame, a local single-particle potential (i.e. a ./V-body non- interacting system) 
that can reproduce exactly the density pi„t of the interacting system. We develop pi n t in the corresponding basis 
{vlnfi * = l-.iV} of one-body orbitals expressed in c.m. frame 

N 

i=l 

We refer the reader to Ref. [T3|, §111. C, for a justification of the introduction of N orbitals in the KS scheme, even 
if only (N — 1) coordinates are sufficient to describe internal properties. (We implicitely supposed that the particles 
are Fermions, but a KS scheme to describe Boson condensates can be set up similarly, equalling all the f\ nt .) 
The KS assumption implies <p\ nt [Pint\ [H|; hence, we can rewrite Ei nt as fl5j 

N - P 2 f 

E lnt [p ln t} = ^(<p'l nt \ — \<p l mt ) + E H [p int }+ E xc [p mt }+ dr v int (r) p int (r), (8) 

i—l 

N 2 

where we have introduced the " non- interacting" kinetic energy X)i=i( < Pintl 2m \f\nt)i ^ ne internal Hartree energy 
Eh [Pint] = \ Jdr dr' p int (r) p int (r') u(r — r') and the exchange-correlations energy functional 

Exc[Pint] = E po tXc[Pint] + EAkin[pint], (9) 

composed of 



E P otxc[Pint] = ^ J drdr ' 7mt(r, r') u(r - r') - E H [pi n t], 



N-l 9 N 



,2 



E Akm [p ln t} = E ^Wint) ~ £>Ll£-|v4t) (10) 

a=l ^ a i=l 

f N p 2 N p 2 

= dr x --- dr N S(R) tp* nt (rx,. . . ,r N ) ^ ^T^( r i' • ■ • > r ^) ~ zJ^ntl^l^"*)- ( u ) 

i=l »=1 

Epotxc contains the quantum exchange and correlations effets linked to the interaction u (7i„t(r, r') is the local 
part of the 2-body internal density matrix defined in Ref. [15j : it is trivially a functional of Pint)- 

In the following, we call "interacting" kinetic energy the kinetic energy of the self-bound system, i.e. 

J dri ■ ■ ■ e?rjv(5(R) ^* nt (ri, . . . , rjv) X^=i ^■V'in^ri, • • • , rjy), and "non-interacting" kinetic energy the kinetic energy 

of the KS system, i.e. Y^i=i( ( Pint\2m\ l Pint)- ^ e see f rom Eq. pT|) that E& kin contains the exchange and "standard" 
correlations 2 , but also the c.m. correlations (due to #(R)), contained in the "interacting" kinetic energy term. It 
is the only term of the functional that contains explicitely the c.m. correlations and represents the main difference 



1 More generally, we can introduce a S(R — a) where a is an arbitrary translation vector, which would lead to perfectly equivalent results. 
We chose a = for simplicity, so that the formalism is formulated in the c.m. frame. 

2 "Standard" correlations means all the correlations except the c.m. correlations in the following. 
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with traditional DFT. Varying Ei nt [pi nt ], Eq. (JS)), with respect to (p\ nt , and imposing normality of the {(p\ nt } leads 
to "internal" KS equations 

/ ti 2 \ 

{ - + U H [pint] + U XC [pint] + V int ) V\nt = WlnV ( 12 ) 

with U H [pint]{r) = 5E H [p m t]/Sp in t(r) and Uxc[Pint}{^) = 5E X c[Pint]/Sp ln t{r), which are local as expected. Equa- 
tions (| 12[) have the same form as the traditional KS equations formulated for non-translationally invariant Hamiltoni- 
ans [l9|, but we have justified their use in the cm. frame for self-bound systems described with translational- invariant 
Hamiltonians, and shown that the functional form of Uxc[Pint] differs by the inclusion of cm. correlations [l5j . 

We see from Eq. (|lip that one has to be cautious with the meaning that is given to the non-interacting kinetic 
energy in mean-field like calculations. Indeed, the non-interacting kinetic energy cannot be considered as a first order 
approximation of the interacting kinetic energy in the general case. The difference is equal to EAkin, that can be 
large when cm. correlations effects are strong, i.e. for small and intermediary size self-bound systems. For large 
self-bound systems E Aki decreases (in relative value), so that the non-interacting and interacting kinetic energies 
become closer. 

The internal DFT formalism has been generalized to time-dependent self-bound systems in Ref. [ill, for instance 
for the description of the collision of two nuclei or laser irradiation, and to multicomponcnt self-bound systems in 
Ref. [l7|, for the description of self-bound systems composed by different kinds of particles (atomic nuclei, mix of 
3 He and 4 He droplets, and molecular systems where the nuclei are treated explicitely) . This last work permits to 
recover the traditional DFT formalism when one kind of particles is much heavier than the others [l7j . underlining 
why traditional DFT is essentially suited for the description of electrons (only) in molecular systems, but not for the 
description of self-bound systems. 

B. A proposed method to obtain a cm. correlations functional. 

We first consider a system of identical Fcrmions (result for bosonic systems will be derived in S jlll F|) . We mentioned 
that the E^ in [p int \ term defined in Eq. (|10[) contains the quantum exchange and "standard" correlations (remind 
footnote [2]) contained in the interacting kinetic term, denoted EAkinXC-, but also the cm. correlations. We thus split 

EAkin as 

EAkin = EAkinXC + E cm , (13) 

where the pure cm. correlations energy E cm is to be defined. Using this definition, the total exchange-" standard" 
correlations energy functional Exci defined by Eq. ([9]), can be rewritten 

rp jpstand i tp 

EjXC = i-XC + ^cm, 

where the so-called "standard" exchange-correlations total energy is defined as Exc 1 = E pot xc + EAkinXC- This 
term is mostly described by the parametrized functionals commonly used for mean-field like calculations of self-bound 
systems (see Refs. 0, Q for a description of functionals used for nuclear systems and Ref. [13] for a description of 
functionals used for Helium droplets systems). But the cm. correlations effects (described by E cm ) are by construction 
mostly not taken into account in those functionals (except through a rcnormalization of the mass in the kinetic energy 
term), which can affect the results, see Ref. [ll|. Our goal is to set up a general well founded form for E cm , which can 
be used for the description of all self-bound systems by simple addition to the common functionals (which, rigorously, 
implies a refitting of those functionals), and is numerically manageable. The idea is to start from the definition (|10p 
of EAkin an d neglect all the exchange and "standard" correlations terms. Then, by definition (fT5|) . we will be left 
with E cm . 

In the following, we propose to search a form for E cm that is functional of the KS orbitals <f\ nt - We follow this 
way because it allows more flexibility, while being fully coherent with DFT. Indeed, the KS orbitals are functionals 
of the internal density, i.e. ip\ nt \pj n t\, as soon as they satisfy KS equations plj , which is the spirit of the Optimized 
Effective Potential (OEP) method [28H30j]. 

C. The commonly used form for the cm. correlations functional. 



We rewrite EAkin[pint), Eq. (ITU1) . in the following perfectly equivalent way 

EAkin [Pint] 



(14) 
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a—1 i—1 



„ * N - 1 2 W 2 

a— 1 " i—1 

« ^ iV 2 p2 ^ 2 

/ dn •• • dr w (V^VintCri, ■ ■ • , rjv)) (£ §: " wv) ^^(n, ■ • ,r ff )) - E^ntl £>Sn*)- 

^ — i — i 



•^/(J(R)V'mt(ri, ■ • ■ , rw) is interpreted as the cm. frame N-body "wavefunction" (remind that ipi nt has the dimension 
of a (N — l)-body wavefunction, see Eqs. (j3|) and (j6|)). This "wavefunction" is obviously not translationally invariant 
(the 5 (TV) hxes the cm. in position space and amounts to move in the cm. frame) and is antisymmetric under the 
exchange of two particles (as ipint is antisymmetric). It is non null only for the that satisfy R = $D i=1 r « = 0, so 
that the {rj} become the cm. frame coordinates. 

Within the internal DFT formalism, the commonly used approximation to treat the cm. correlations can be 
recovered by supposing that the KS Slater Determinant, denoted ip aux , is a good hrst order approximation of the cm. 
frame N-body "wavefunction": 



VWfe(ri,..,r w ) « ^(n,...,^), (15) 
where (P are the possible permutations of the coordinates and p the number of transpositions of P) 

r(n, ...,r N ) = ^= E(-i) p nf =1 ^«(r 4 )- 

V ivl p 

Inserting this approximation in (fH)) and following the method described in ijll Bl (the " standard" correlations are by 
construction neglected and the exchange terms cancel naturally), we obtain 

E Akin ~> EcnilWint}} 
I I P 2 I 

V \2mNV 



h 2 N r r 2 N r r 

= -^E / dr Ct(r)Ar^ n( (r) - 7^ E / *^WV,4(')/ rfr'^ 4 (r')V r ^L t (r'). (16) 

i—1 i.j — l 

We recover the commonly used form for the cm. correlations functional (note that, in practise, the second term of 
the last line of Eq. (fT6|) is often neglected to reduce the numerical cost) [ll| . 

The internal DFT formalism permits us to shed a new light on the validity of the approximation ([TBI) . It holds if 
and only if the approximation (IT51) holds at least to first order. But this cannot be the case in general because ip aux is 
far form being null when J2i=i r i ^ 0- Moreover, ip aux contains a cm. vibration, typical of Slater determinants (i.e. 
^a«x|pn|^«x) q for n > 2 ) [l|,[|, whereas s/S(Tl)tp int does not contain such a vibration (i.e. (^™*|P n |V> int ) = 0, 
Vn). Thus, we cannot expect for a systematic satisfying improvement with this form [11] . 

In the next section, we propose a better founded form for the cm. correlations energy functional, where the cm. 
correlations (the 8(TV) term) appear explicitely. 



III. A GENERAL NEW FORM FOR A LOCAL CM. CORRELATIONS POTENTIAL. 



A. The idea. 



We adopt a different point of view than the one of mi CI We start from E &kin[pint] written as in Eq. (fTTj) (instead 
of (HI])) and pose (instead of (|T5]|) 

1 p int (r 1 ,...,r N ) -+ —t— r ™( ri ,...,r N ). (17) 

The reasons of the introduction of the 1-body "wave- function" r aux (TV) (that implicitcly depends on the number of 
particles N) are the fallowings: 
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• ipint has the dimension of a (N — l)-body wavefunction, whereas the KS Slater Determinant ip aux (ri, . . . , rjy) 
has the dimension of a TV-body wavefunction. Dividing ip aux by Y aux permits to recover the dimension of a 
(N — l)-body wavefunction while preserving antisymmetry. 

• The KS Slater determinant %\j aux contains a cm. vibration, as already mentioned. This is not a problem from 
the KS point of view, because all is wanted is that ip aux reproduces the correct pi n t , the KS Slater determinant 
representing nothing more than an auxiliary. Hovever, if we wan't to use tp aux to approximate if>i n t, the cm. 
vibration of ip aux has to be "substracted" , because ipi n t should not contain such a vibration. r aux represents 
the proposed way to "substract" the cm. vibration contained in ip aux . 

In the particular case of the harmonic oscillator (i.e. when the interaction u is parabolic), we can al- 
ways achieve the separation vp aux (r 1: . . . , rjv) = r aii:E (R) x F(£ 1; . . . , £,n-i) EH- Thus, the use of the form (fTTj) 
permits to "substract" exactly and directly all the cm. vibration contained in ip aux . 

In the general case, ip aux (ri, . . . , rjy) cannot be separated into r otla:: (R) x . . . ,£jv_x). However, all 

we wan't is that the form of T aux would allow sufficient flexibility to lead to a satisfying result for both E cm 
and pint (then, even if the " substraction" is not "direct" because there is no separation of the cm. motion in 
ip aux , it is "indirect" because it leads to the correct final result). 

• As we will see, the final result has a clear Physical meaning, which shows the pertinence of our proposition. 

• The 5(R) contained in (fTTj) breaks translational invariance, so that even if rau l^ R s ) ip aux (ri 7 . . . , rjv) is not 
translationally invariant (contrary to ipint), the final result won't be affected if a good form for Y aux is found. 

Necessary properties that T aux should satisfy will be discussed in mil Dl 



B. The result. 



We first note that, when inserted in (fTTj). the approximation (fT7|) does not necessarilly lead to a real E Akin 3 . To 
overcome this problem (which would affect time-dependent results), we insert (I17|) in E Akin and keep only the real 
part. We then obtain an approximation of the exact E Akin where the "standard" correlations have been neglected 
by construction. As discussed in £111 B( it remains to neglect the exchange terms in the result to obtain E cm . The 
calculation is detailled in Appendix [A] The final result is 

E Akin -> EcrniWint}} = 

fc2 ^ r. AC 

E J rfr vfoW^vLM >< ( jrss(5jj5 / dr ' x hiM{<p%t' l }}(* + O - 



i=l 

h 2 



2mN 

+ Purelmaginary[{(pi nt }], (18) 
where the functional Purelmaginary is defined as 

fc2 N r i r 

Purelmaginary[{rf nt }] = / dr ^WA r ^„ ( (r) x J dr' |^(r')| 2 x /i, ¥ i[{^}](r + r'). 

(19) 

It counter acts the pure imaginary part of the first line of ([T8|l and becomes null in the real (stationary) case. We 
keep it for generality (time-dependent case). 

The "two particles cm. correlations functional" is defined as (D is the dimension in which the calculations are 
done: D=l,2 or 3) 

hiHirft: ' l W) = n d I n*L dr,s( J2 r k + f)n* =1 \<f 



3 This is fundamentally due to the fact that, in the general case, the form H17I I cannot be rewritten as a function of the only. 
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fi i^ti couples the particles i and I together in (|18p . The meaning and properties of this functional will be detailed in 
The potentials U l cm corresponding to E cm are defined by 4 



~£{ Ar ^ (r) x i\r^W / *^W'A>^^' m M T + ^ - l ) 
i N r 

+ | r an, (0) |2 ^( r )E 7 dr ' ^ t (r')A r ^L t (r') x / W i[{^}](r + rO} 

2mN K K^^ AR r^R)lR=o) X ^ (r) / dr 'l^( r ')! 2 X funMW-^ m }}(r + r') 
5_ 



+ c ;* f ^ PureIr naginary[{ipf nt }}, (21) 



where the last line is obviously null in the real (stationary) case. Note that the potentials U l cm are not common for 
all states. This is due to the fact that E cm [{ip^ nt }] is orbital dependent, which necessits extra measures to recover a 
common potential or, equivalently, to preserve orthonormalization. A way to overcome this problem is to use the OEP 
method, which permits us to find the potential common to all states that reproduces the most accurately the effect 
of the U [ cm potentials. We refer the reader to Refs. (28l - [30j for the exhaustive equations. As the full OEP result is 
very costly numerically, it is often simplified. Popular is, e.g., the Kriegcr-Li-Iafrate (KLI) approach and, in a further 
step of simplification, the Slater approximation [3lT [32j . As our goal is to find a numerically cheap form for the local 
cm. potential, we detail thereafter only the Slater approximation. It gives: 



' l^(r)| 2 
^ /w(r) 



U!L a \r) = £ Wte&LuU')- (22) 

i=i 



C. Properties of fij^ti. 

The definition ([2H)l of the two particles cm. correlations functional fi^i shows that: 

1. it is real and has the dimension of a density, 

2. it is normalized to N D , i.e. J dvfi^iii) = N D , 

3. lim^ioo /i,{^i(f ) = 0, 

4. it is a " multiconvolution" of all single densities, unless these associated to orbitals i and 

The first three points permit to explicit the physical meaning of jp? fi,i^i(r + r'): it is the probability that particle 
I =/= i has position r', given particle i has position r. Indeed, because of the cm. correlations, the positions of 
those particles arc not independent and are coupled by fi.i-a- As a consequence, this coupling appears in the cm. 
correlations energy (fT8|) and potentials (|2~Tj) : because of the cm. correlations, every single orbital (p\ nt is coupled to 
every single orbitals ip 1 ^. through fij^i. 

To understand better this coupling, note that /u^i can be rewritten 

hiMWt^W) = 2^(r) if N = 2, 

3>**''(-f)| a if ^ = 3, 



4 A usefull and straightforward relation to achieve the variation is (<5; ; is the Kronecker symbol): 
j^fK IfZf (r')l 2 x/ Wi [{^„1 i ' m }](r' + r") = (l-<5 M ) vLtM fuMto^K* + r"). 
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N 

N D n N ^ dr,U N 3 ^ |^„ t (r,)| 2 |^ t (- E r *- f )! 2 ^ ( 23 ) 



Constant, for very large N (limit of a Fermi Gaz). 

We see that, in the 2 particles case, f\$ is proportional to a the steep delta function. Thus, if particle 1 has position 
r, particle 2 will have position -r, so that the cm. remains stuck to R = 0. In the 3 particles case, f%i^% has a larger 
width, because the introduction of a third particle allows more freedom to the motion of the 2 other particles, while 
preserving R = 0. As /u#i is a multiconvolution of all the single densities unless the ones associated to orbitals i 
and I, it should fastly become very little dependent of i and I as N grows. For N > 4, the width of will increase 
as N grows, because of the multiconvolution form of fn^i- Indeed, a larger number of particles allows more liberty 
to the motion of 2 of these particles, while preserving R = 0. For very large N, the system tends to a Fermi gaz, so 
that fij^i tends to become constant and dclocalizcd in the whole space, i.e. the motions the particles tend to become 
independent. The cm. correlations can then be neglected, as expected. 

Practically speaking, we see that the numerical cost of the whole scheme lies in the calculation of /u^i for N > 4, 
i.e. the calculation of the multiconvolution of Eq. (|2"3"|) . At first sight, it seems to be discrimatory for large N. 
But a mathematical property of the convolutions under Fourier transforms makes it manageable. In Appendix [B] is 
demonstrated the so called "multiconvolution theorem". Its direct application to for A > 4 gives (T denotes 

the Fourier transform as defined in Appendix [0 Eq. (|BT])) 

This permits to fasten drastically the numerical calculation of fi^i, which becomes manageable even for large systems. 
Indeed, once every T[|(/?i nt | 2 ] are calculated, fi^i is given by the inverse Fourier transform of their direct product, 
so that the numerical cost of fijjti equals the numerical cost of (N + 1) Fast Fourier Transforms when N > 4. 



D. Properties of T aux . 

Remind that Y aux implicitcly depends on N . We see from (fT8"|) that Y aux should be two times derivable and that 
only the values of Y aux and of its second derivative in R = are of interest to reproduce E cm . We now characterize 
those values. 

First, we note that |r au:r (0)| 2 is imposed by the normalization condition on the approximation (fTT)) we used for 
ipint 5 : 



1 = (iptnt\ipint) = J dri ■ ■ ■ dr N 5{B.)\ip int (r 1 , . . . ,rjv)| 

=> |r— (0)| 2 = fdr dv' |^„ t (r)| 2 |^(r')| 2 /u^[{^1 M }](r + r')- (24) 



Numerically speaking, this condition will be satisfied self-consistently, starting from a reasonnable initial value for 
|pau:c^Q-j|2 an( j resca i m g ft a ^ every loop so that it satisfies (f2"4"]l . 



It remains to characterize the value of rat ,i» Ar 



i a_ i 



strated in Appendix [Cl 



. To that aim, we use the following results demon- 

R=0 



lim \T aux (0)\ 2 -> +cx), (25) 

Af->+oo 



lim — x — 



n^+oo N r aux *(0) r aux (R) 



x \r aux (0)\ z -> 0. (26) 

R=0 



The most simple and natural [J| form for Y aux that meets those criteria (and is "exact" in the case where the 
interaction u is parabolic) is a Gaussian 

i=l 



5 The second line of Eq. i 124[ ) is obtained introducing the form H17H for 4'int in the first line of Eq. i l'24[ ) and neglecting the exchange terms. 
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where Ri are the coordinates of R in D dimensions and K(N) is a two times derivable function of N. With this form 

|r— ( )| 2 = (Wh D/ \ (28) 



7T 



1 -A R1 1 



paux*(Q) paws fj^) 



/ x l-D/2 

= tt d/2 x D x X(iV) . (29) 

R=0 V / 



This permits us to define completely the cm. correlation energy (fT8"| and potential (f2~Tj) . without introduction of 
any free parameter. Indeed, as previously mentionned, |r aua: (0)| 2 is defined by Eq. (f2"4"|) and should be rescaled 
at every numerical loop so that it satisfies this relation. This defines K{N) at each loop by (|28[) . which defines 

r a„i»( ) AR ra „^ R j by (f2"9")l , and we are done. 

With this method, there is no need to define analytically K as a function of TV; K(N) is obtained numerically for 
every given N as indicated. It would nevertheless be interesting to obtain an approximate analytical form for K(N), 
at least to start the numerical iterations with a pertinent initial condition. To that aim, we notice that the conditions 
(f25j) and (f26|) . together with the equalities (|28| and (|29|) . imply the following conditions on K(N) 

lim ( K(N) ) -> +oo, 

lim — A'(iV) 0. (30) 

jV-f+oo N 

A straightforward form for K(N) that satisfies those two constraints is 

K{N) = Ax N a , where < a < 1. (31) 
The corresponding energy associated to T aux is then 

ti 1 h 2 n h 2 n 
E rau * = -__(r aua: |A R |r ara ) = ———K(N) = - -A x /V 1 " 1 . 

Ey°-w is then proportional to N 11 ^ 1 . Even if has, strictly speaking, no Physical meaning, is is reasonnable 

to assume that its variation according to N should approximately be proportional to the variation of the energy 
associated to the cm. vibration obtained in mean-field like calculation, see Ref. For the nuclear case, the cm. 
vibration energy evaluated for 3D harmonic oscillator states is proportional to TV -1 / 3 ; this variation can be reproduced 
with a = 2/3. The cm. correlation energy evaluated with a a posteriori fit with 3D mean-field like calculations is 
proportional to N~ - 2 ; this variation can be reproduced with a sa 0.8 [ll|. This is fully coherent with the fact that 
< a < 1, as deduced by our previous considerations, which reinforces the pertinence of the Gaussian form, that we 
will use in the following numerical calculations. However, the search for other forms, i.e. with other variations around 
R = 0, should be continued to obtain the most precise description of self-bound systems in fully realistic calculations. 



E. An explicit density functional for the cm correlations energy. 

The functional proposed in ^IIIBl is well suited for stationary calculations, but not for time-dependent ones, because 
the Slater approximation, as well as the KLI approximation, do not permit to preserve energy conservation (as those 
approximations are not perfectly variational }33l|). Only the full time-dependent OEP result [34| will achieve energy 
conservation, but at the price of a much larger numerical cost. It thus would be interesting to find a functional that 
overcomes those pathologies in the time-dependent case. 

In this section, we propose a further step of approximation that will allow to obtain an explicit functional of pi n t 
for the cm. correlations energy. We will do the Local Density Approximation (LDA) on the result of ^IIIBl which 
consists in doing the approximation that the system is locally homogeneous [U |35| . Despite of its simplicity, this 
approximation has proven to be very satisfying for the description of a wide range of systems and not only large 
ones dlUl. To do the LDA, we first do the replacement 

l^ntMl 2 -> ^Pint(r) (32) 

everywhere the single density terms \tp\ nt \ 2 appear into E cm , Eq. (ITS)) . In the obtained functional, the only remaining 
term that is not an explicit functional of pint is X^=i ^* t (r)A r <p| ni (r). As we consider a system composed of Fermions, 
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we do the Thomas- Fermi approximation [21[ , i.e. the replacement 

N „ 

5>» t (r)A r ^ nt (r) --Cp^ir), (33) 

i— 1 

where C = (^-) 2 ^ 3 and 7 is the degeneracy 6 . We obtain as a final result the cm. correlations energy written as an 
explicit functional of Pint'- 



E^ A \pint\ = ^Jdr §<7p& 8 (r) ( |rJ (0) |2 / x M/^Kr + r') - l) 

3f?ef . ; Ar 5—— \x f dr —p, n t(r) [ dr' — p int (r') x h\pi n t\{r + r'), 

\Y avx *(0) T aux (R) R=0/ / jyruiLy 1 1 jyrwM\ y jA\s-mi,i\ n 



(34) 



2to7V 

where the "two particles average cm. correlations functional" is defined as 



1 f N ~ 2 

fe=l 



The corresponding unique cm. correlations potential is given by 7 



W = ^ MW = £{^ (r) X (^4h)F / ^ + O - l) (36) 

h 2 



5ie( ^— -A R \— )x [ dr'—p int (r') x f 2 \pi n t}(r + r'). 



2mN 

fi has the following properties: 

1 . it is real and has the dimension of a density, 

2. it is normalized to N D , i.e. / df/ 2 (f) = iV D , 

3. lim f ^ ±oc / 2 (f) = 0, 

4. it is a " multiconvolution" of (N — 2) /jj„t. 

The first three properties and the considerations of 3111 CI permit to interpret 5^7/2 (r + r') as the average probability 
that one particle has position r', given another particle has position r. / 2 can be rewritten more explicitely in the 
following way: 

f 2 \pint}(v) = 2 D S(r) if N = 2, 

3^ mt (-r) if N = 3, 

JV-3 

D 



^—^ f drx . ..dr N _ 3 p int {r{) x ■■■ x p int {r N - 3 ) p int (- ^ r fe - f) if N > 4, 

^ k=l 



Constant, for very large N (limit of a Fermi Gaz). (37) 



6 More precisely, if we take into account the dimension D in which the calculations are done, we have 

££1 ^» t (r)A r ^ nt (r) -JB-C D p*f (t), where C 3 = (^) 2 / 3 , C 2 = f and & = 2=i. 

7 A uscfull and easy to demonstrate relation to achieve the variation is the following: 
^if^y/*' jrPint(r')xf 2 [p int )(r' + v") = V \ nt (r) »=k / 2 [ft„t](r + r"). 



11 



Still, we see that the numerical cost lies in the calculation of f2[pint] for TV > 4. To reduce this cost, we use the 
" multiconvolution theorem" demonstrated in Appendix [Bj Its direct application to f2[Pint] for TV > 4 gives (using the 
definition Eq. (|Blj) for the Fourier transform T) 

fa[p in t}(r) = N D j^ x T- 1 [{T[ Pmt ]) N - 2 \ {-¥). 

This permits to fasten drastically the numerical calculation of fa, which becomes manageable even for very large 
systems. Indeed, one simply calculates T[pi n t], puts it to power (TV — 2) and calculates its inverse Fourier transform. 
Thus, the numerical cost of the calculation of fa is equal to two Fast Fourier Transforms for all TV > 4. 

Moreover, this scheme is perfectly variational, contrary to the one of £ )III Bl as already mentioned, and thus will 
be suitable for stationary calculations as well as for time-dependent ones (i.e. it will achieve energy conservation if 
time- independent Vi n t is used). 

F. And for Bosons? 

The cm. correlations energy functional for Bosons condensates is obtained simply by doing f\ nt — > <J>int, thus 
Pint = N\ip mt \ 2 , in (HHJl. Posing (p int (r) = y/p int {r)/N } we obtain 

(38) 



E cm [ptnt] = J dr y/ pint(r)A ry / p int (r)x [ p aux ^ 2 J dr' — p int (r') x fa[p int ](r + r') - 1 

^ e ( T auL( 0) A R raM ! (R) | R=0 ) x J dr ^P™t(r) Jdr' ^ Pmt (r') x fa[ Pint \(r + r'), 



2mN 



where the two particles average cm. correlations functional is defined as in Eq. (|35|) and (|37|) , and the corresponding 
cm. correlations potential is given by 

5 fptt{rt =Ucm[Pmt](v) = "t{7^ ArV ^K^0F/ dr '^ mt(r ' )/2 ^ 

+ |T°^(0)p / dr ' VPmt(r)A r ^p mt (r) x — ^fa[p int ](r + r')} 
H^^ AR f^R)L=o) X Jdr'^p mt (r')xfa[p mt }(r + r'). (39) 



2mTV 

This potential is common to all states, so that no extra measures to preserve orthonormalization are to be taken, and 
it is strictly variational, so that it is usable in the time-dependent case. 

IV. FIRST NUMERICAL RESULTS. 

In the following, we will consider ID calculations, which permit to understand better some features of the internal 
DFT formalism and to include more easilly various particle-particle interactions. 

A. Model system composed of two different particles with a strong interaction. 

1. The model and the benchmark. 

We consider a stationary multicomponent self-bound system composed of two different particles, to maximize the 
cm. correlations effects. We suppose that the two particles have the same mass m and are coupled by a strong 
interaction, which models for instance a proton and a neutron. The first particle has laboratory coordinate A 1 * 1 and 
the second has laboratory coordinate A 2 \ The reduced mass is p, = m/2 and the Jacobi coordinate £ = A 1 ^ — A 2 \ 
We suppose that the interaction between the two particles is parabolic (harmonic oscillator) , so that the Hamiltonian 
is 

„(1) 2 r,(l) 2 1 
Tj P i P i 1 2/(1) (2)\2 

H = — h — h -muj (r K ' ~ r x ') 

2m 2m 4 ' 
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The exact solution can be written in an analytical form (ipint should not be anti-symmetrized because we deal with 
two different particles) 

The corresponding energy is Ei nt = \hx> and the cm. frame one-body densities for each kind of particle are fl7| 
(R = ( r W +r (2))/2): 

plStW = /"«irWdr«<y(ii)|V<nt(rW-r( 2 ))| a *(r-(rW-ii)) 



= 2|^ ( 2r)| 2 = ^exp{-^r 2 }. (42) 

This is our benchmark. 

It can be shown analytically (using a harmonic oscillator basis) that the Hartree (H) solution (remind that there is 
no exchange because the two particles are different) leads to E int = and /oj„ t (r) = J exp(~-^=|r 2 ). Thus, 

the H energy is 2/^/2 (sa 1.4) times more important than the one of the benchmark and the density is 1.7 times more 
spread. The H solution is much more delocalized than the benchmark because the cm. correlations are neglected 8 . 



2. The internal DFT exact functional. 

Applying the multicomponent internal DFT formalism developped in Ref. jjjj (whose equations have a relatively 
similar form than the " one kind of particles" internal DFT ones recalled in mi A[) we can rewrite the internal energy 

as a functional of p^ t and p\^ t (we have choosen v int = 0) 

F n {2) ] - (J 1] \ p2 \in {1) \ + (in {2) \ p2 \in (2) \ + F (12) \n (1) n (2) ] 4- F (12) b, (1) n {2) ] 4- F {12) \n (1) o {2) ] (M\ 

^intlPinV Pint! ~ VPint I ^ I Vint I + l^int I jj^ I Pint ) + ^ H [Pint' Pint! + ^potclfiinV Pint! + ^AkiniPinV Pi-nth V i6 ) 

where 

E ( H ] [pflp? n \] = j drdr'p^ t {r)pf^)\ m u 2 (r-r>) 2 (44) 

F {12) r (1) (2)i _ / , iT 2 , , , (1), P 2 , (IK r (2) | P 2 , (2)x 
^AkiniPint'Pinti ~ \Vint\-^\Vint) ~ Wintly^Wint) Wtntl ^ Wint )■ 

(12) (12) 

E y H 7> is the H energy, E potc is the "standard" correlations energy linked to the particle-particle interaction and 

(12) 

^Akin i s t ne energy associated to the correlations contained in the interacting kinetic energy; it is the only term that 
contains explicitely the cm. correlations. Using the results of flIV A 11 we can reformulate those functionals as explicit 
functionals of p\ n \ and p^ t . We have [13] 

-Yint\rS) = f dr^dr^5(R)\^ mt (r^ - r^)\ 2 6(r - (r« - R))S(r' - (r^ - R)) 



The cm. correlations tend to localize the densities compared to the independent particle approximation, which can be understood as 
follows: if cm. correlations are taken into account, when one particle moves in a direction where the potential well is higher, the other 
one will have to move in the opposite direction where the potential well is also higher; thus, the first particle will feel the repulsion 
present in the independent particle approximation, but also the repulsion that feels the second particle through the cm. correlations, 
which will localize its density. 
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M mt {2r)\ 2 5{r + /) = p^ t (r)S(r + r') = \(p^ t {r) + p<&(r))s(r + r'). (45) 



Thus (we pose ip, 



(0 



Pint ) 



p (12) r (1) (2)n . 1 ( n W( r \ , n (2)^ ffl , ,2_2 F ( 12 )rJ!) J 2 )l 

"f^ + l / dr (^nt( r ) +Pi 2 l t( r )) mw2r2 - 



K (12) « (2) 1 
^Akin [Hinti Pint) 



(46) 
(47) 



3. TTie cm. correlations functional. 

1/4 



We now apply the idea presented in ^IIII with T aux (R) = cxp j — -|-i? 2 j. We obtain 

^ = ~ / *[eMA^ t W(2 J|l^l(-r)| 2 - l) + eWA r ^ t W(2^|^(-r)| 2 - 1 



y ^i^wn^c-oi 2 , 

and for the local cm. correlations potentials (I = 1,2) 



^« ^ -^£( 2 ^ l ^' )( -*- )|J - 1 ) A '- v! "' (r) 



2m 



V V # (- r ) A ^mt (-H +vX7r|^„^ ; (-r)| 



(48) 



(49) 



^. Numerical results. 

We use a unit system where K = m = 1 and we choose w = 1. Table U and Fig. [TJ give numerical results for the 
following formalisms: 

• Benchmark (described in ^IV A 1|) . 

• e£ 2 \ called "H only", 



• E 



• E 



• E 



• E 



(12) 
H 

(12) 
H 

(12) 
H 

(12) 
H 



(12) 

E y c ' , called " H + standard correlations" 



E, 



(12) 



- < 



2iiiN 



>, called "H + standard correlations + standard cm. correction", 



E. 



(12) 



41, called "exact internal DFT", 

E cm , called "internal DFT with cm. correlations functional"; we obtain K w 3.9 with the 
method described in ^III Dl 

First of all, we sec from table U and Fig. [I] that "exact internal DFT" reproduces perfectly the total energy and the 
densities of the benchmark, so that the non-interacting v-representability [TH, [2l| is perfectly achieved. This is not a 
surprise because the non-interacting v-representability is trivially achieved when one deals with only one particle of 



E. 



(12) 



each kind (because it is always possible to reach tp 



(0 _ JM ) 

~ V Hint J ■ 



Note that, even if in the general case the functional E^^-, + E^^ n is universal |15|, |16|] , the forms I I46II and J 471) are limited to the two 
different particles case because exchange effects are not taken into account; they thus cannot be used to describe a system composed 
by and arbitrary number of particles of each kind. The universal functional, applicable to an arbitrary number of particles, is more 
involved but should permit to recover (146 P and II47I I in the limit of a system composed by two different particles. 
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Formalism 


Non-interacting kin. energy 


K (12) 

H 




- < 2^V > 01 E Akl 01 E c,n 


Total energy 


H only 


0.353 


0.353 








0.71 


H + stand, corr. 


0.5 


0.25 


0.25 





1.00 


H + stand, corr. + stand, cm. correct. 


0.706 


0.177 


0.177 


-0.353 


0.71 


Exact internal DFT 


1.000 


0.125 


0.125 


-0.750 


0.50 


Internal DFT with cm. corr. ft 


1.225 


0.120 


0.120 


-0.918 


0.55 



TABLE I: The energies of the various formalisms (in units where h = m = 1; benchmark: total energy = 0.50 and interacting 
kinetic energy (ip in t ||-| ip int ) = 0.25). 




FIG. 1: (Coloured version online.) The internal densities p^t °f the various formalisms (x-axis: position in units where 
H = m = 1). 



From table HI we see that the non-interacting kinetic energies of the various mean-field like calculations cannot 
be compared to the interacting kinetic energy In particular, there is a factor 4 between the "exact internal DFT" 
non-interacting kinetic energy (equal to 1) and the interacting kinetic energy (equal to 0.250). Indeed, as already 
discussed in £111 Al it is the "non-interacting kinetic energy + E Akin " that is comparable to the interacting kinetic 
energy. Exact internal DFT then reaches perfectly the benchmark: 1.000 — 0.750 = 0.250. Internal DFT with E cm 
gives 1.225 — 0.928 = 0.307, which reproduces fairly well the benchmark considering that the "standard" correlations 

(12) 

part of E Akin has been neglected. The result with standard cm. correction gives 0.706 — 0.353 = 0.353, which is 
worse. 

From the point of view of the total energy, "H only" is far away from the benchmark, as well as "H + stand, corr." 
and "H + stand, corr. + stand, cm. correct.". "Internal DFT with cm. corr. ft." is the closest to the benchmark. 

From the point of view of the densities, Fig. [1] shows that "internal DFT with cm. corr. ft." is very close from the 
benchmark and represents a great improvement compared to the other results. Nevertheless, we see some differences 

(12) 

that should mostly be due to the fact that the "standard" correlations part of E Akin is not taken into account. In 
realistic mean- field like calculations, however, this "problem" should not appear because those "standard" correlations 
are mostly taken into account in the parametrized functionals that are commonly used. 

B. Model system of identical Bosons (smooth interaction). 



We now study the results of internal DFT on another ID system, with a "smooth" interaction. We consider a 
system composed of N identical Bosons of mass m and positions {r{}, without spin, in a condensate state and with 
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an attractive 2-body interaction of the form (e > 0) 

u(r — r ) 



yf{r-r') 2 + e 



(50) 



The more e is large, the more the potential is smooth. This allows to model, in some way, some features of 4 Hc 
droplets. 

The internal DFT energy functional is given by (pint is the 1-body orbital describing the bosons with pi nt = N\(fi nt \ 2 
and we consider Vi„t — 0) 



„2 i 

E in t\Pint] = N(if int \ — \ip mt ) + E H [p mt ] x (1 - — ) + Eg and [p int ] + E cm [p int ], 
Ira N 



(51) 



where (1 — jT)Ejj[pint] represents the H energy where the self-interaction has been substracted and E^ and [pi nt ] is the 



"standard" correlations energy defined as 

If 1 

Ec tand [pint] = - / drdr'-y int (r,r')u(r-r') - E H [p mt ] x (1 - — ) + E Akin [p int ] - E cm [p mt }. 

The cm. correlations energy is defined as in mil F[ where r aux is defined as in mil Dl 

dr y/ Pint{r)A ry /p int (r) x ^^N) j dr ' ]^«*( r ')/2[Pint](r + r') - l) 



E cm [Pint] 



2m 

h 2 



W V^K(N) I dr \^ mt (r)\ 2 I dr' ^p mt (r') x f 2 \p int ]{r + r'), 



and ji is defined by Eq. (|37|) . The internal KS equation is 

/ h 2 1 \ 

( - —A + U H [pint] x (1 - — ) + Ug and [p int ] + U cm [p in t})(p ln t = epint, 

where Ug and = SE s c tand /5p mt and 

U cm {p lnt ]{r) = -^-\ —^—— A r \/p m t{r) x (</^7^r / dr' ^-p int {r')f 2 [pint\(r + r') - 1 
2m l y / p l nt(r) V y K(N) J N 

KjN) j dr ' ^ P™t{r')A r yJ Pmt (r') x — ^±f 2 [p int ]( r + r')} 



-^V^V> f dr'^p mt {v') x f 2 [p mt }(r + r'). 



(52) 



(53) 



(54) 



Formalism 


Non-interacting kin. energy 


E H x(l- i) 


" < > or E om. 


Total energy 


H only 


0.133 


-0.626 





-0.49 


H + stand, cm. correct. 


0.260 


-0.712 


-0.065 


-0.52 


Internal DFT with cm. corr. ft 


0.535 


-0.776 


-0.418 


-0.66 



TABLE II: The energies of the various formalisms in the N = 2 case (in units where h — m = 1; benchmark: total energy 
= —0.59 and interacting kinetic energy (ipintlj^lif'int) = 0.12). 

For the N = 2 case, we can compute a benchmark. Indeed, by use of Jacobi coordinates, the internal Hamiltonian 
can be rewritten H. mt = ^ ~ ~/p = i where // = m/2 is the reduced mass. It is then possible to calculate numerically 

the exact many body solution ground state ipint and the cm. frame one-body density p int (r) = 4|-0.;„ t (2r)| 2 . 

In the following, we neg lect Eg and because we have no simple way to evaluate it as a functional of pint and it 
should be a correction in the general case (even if not always fully negligible). The next results are given for: 
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0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


-4-3-2-101234 

FIG. 2: (Coloured version online.) The internal density pint/2 of the various formalisms in the N — 2 bosonic case (x-axis: 
position in units where h = m — 1). 



N = 


Non-interacting kin. energy 


Total energy 


2 


0.133 


-0.491 


3 


0.385 


-1.75 


4 


0.717 


-3.79 


5 


1.138 


-6.66 


6 


1.632 


-10.37 



TABLE III: The "H only" energies for various N (in units where h = m = 1). 

• E H X (1 - jj), called n H only", 

• Eh x (1 — tt) — < o P at >, called "H + standard cm. correction", 

n V N ' 2mN ' ' 

• Eh x (1 — jf) + E cm , called "internal DFT with cm. correlations functional". 

Table Eland Fig. [2] show the energies and densities in the N = 2 case. We see that "internal DFT with cm. corr. 
ft." reproduces fairly well the benchmark, at least much better than "H only" and "H + stand, cm. correct" (we 



N = 


Non-interacting kin. energy 




Total energy 


K 


Interacting kin. energy 


2 


0.535 


-0.418 


-0.66 


1.94 


0.117 


3 


0.463 


-0.185 


-1.90 


1.74 


0.278 


4 


0.702 


-0.196 


-3.97 


2.74 


0.507 


5 


1.014 


-0.217 


-6.84 


4.01 


0.799 


6 


1.390 


-0.239 


-10.56 


4.54 


1.151 



TABLE IV: The "internal DFT with cm. corr. ft." energies for various N (in units where h — m — 1). 
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0.7 



0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


-4-3-2-101234 

FIG. 3: (Coloured version online.) The internal density pmt/6 of the various formalisms in the N — 6 bosonic case (x-axis: 
position in units where h = m — 1). 



1.8 



1.6 
1.4 
1.2 



0.8 
0.6 
0.4 
0.2 


-4-3-2-101234 

FIG. 4: (Coloured version online.) /2 for N — 3 to 6 (N — 2 is a delta function; x-axis: position in units where h = m = 1. 



nevertheless see some differences, which should mostly be due to the fact that E^ and has been neglected). Again, we 
see that the internal DFT non-interacting kinetic energy cannot be compared to the interacting kinetic energy. It is 
the internal DFT "non-interacting kinetic energy + E cm " (0.535 — 0.418 = 0.117) that is comparable to the interacting 
kinetic energy (0.12), as underlined in £)II Al The result with standard cm. correction gives 0.260 — 0.065 = 0.195, 
which is worse. 

Tables Hill and HVl show that "H only" total energy becomes closer to internal DFT total energy as N grows. 
Moreover, table IIVI shows that the part of the cm. correlations energy in the total internal DFT energy decreases 
as N grows (63% for N — 2; 2% for N — 6). As a consequence, and even is shell effets play a role, the internal 
DFT non-interacting kinetic energy tends to become closer to the interacting kinetic energy as N grows (factor 4.6 
for N = 2; factor 1.2 for N = 6). The "H only" and internal DFT densities also becomes closer as N grows, see Fig. [3] 
for N = 5, as awaited. 

Finally, we see from table HVl that K grows as N grows and from Fig. [3] that the maximum value of also grows 
as N grows (but the maximum value diminishes), confirming the reasonings of ^IIIDI and Appendix [C] 



V. CONCLUSION. 

Internal DFT provides an existence theorem for a cm. correlations energy functional associated to a local potential. 



H only 

Internal DFT with cm. cor ft 
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In this article, we have constructed such a functional, without including any free parameter. The use of this functional 
is justified by a strong formal background and variants suitable for fermionic as well as bosonic systems have been 
proposed. The resulting scheme is numerically manageable and represents a well founded alternative to projection 
techniques to treat the cm. correlations. It can directly be added to actual energy functionals (but a refuting of them 
then should be necessary). Moreover, it permits to recover clearly and precisely the value of the interacting kinetic 
energy and represents a manageable way to include the cm. correlations in time-dependent calculations of self-bound 
systems. 

We have presented first convincing numerical results on ID model systems. They amongst other show that the 
developped functional represents a great improvement compared to the "standard cm. correction" commonly used 



in nuclear Physics (of the form 



< 



2mN 



>), especially from the point of view of the energies. The next step will 



be to include the proposed cm. correlations functional in realistic 3D calculations, for instance in mean-field like 
calculations of nuclei with Skyrme interaction @, Q . As the " standard" correlations are mostly taken into account in 
the commonly used functionals, the 3D results should be even more convincing than the ID ones. 

Finally, even if the proposed Gaussian form for r aux (R) has been proved to be pertinent, it could be interesting 
to search for other forms, i.e. with other variations around R = 0, to obtain the most precise description of atomic 
nuclei, Helium droplets or molecular systems with explicit treatment of the nuclei. 
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Appendix A: Details of the calculation that leads to E C m- 
We evaluate / dri . . . drjv^(R)V'int( r ij ■ • ■ > t n) J2iLi J^rV'rot( r i) • • • , ?n) using the approximation (|T7|) for Vin* an d 



obtain 



10 



f dri . . . dr N S(R) __-^""*(r lj . . . ,r N ) V L_^° u *( ri) ...,r N ) 

N 

E(-!) P+P 'E / *< *f ,p, [{*£*}]('«) >< ^(*o^2%) 

n n, i=1 J 

i N r 

n n / a i J 



(Al) 



V 1 1 

2m |r aua: (o)| 2 aT 
K 2 1 1 1 

'2^ ~fp r aux*(Q) R r ara (R) 



J^E{/*«K*<£?>K') x ^ t (r)A r ^„ t (r) 

p P'^p ^ 

R=0 ^ E { / dr }]W x K lt (r)| 2 

?'— 1 



1 A 1 

"Ar- 



2mA r au ^* (0) r aM:r (R) 



P,P'=£P 



Note that A r 



C (R) 



i 



, which is obtained by a simple variable change and is usefull for the calculation. 
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where we have defined (D = 1, 2 or 3 is the dimension in which the calculation is done) 



Fr p 'iwLw) = n d 



k^i 

and its diagonal part 

JV 



r N 

/ nf =1 dr i 5(^r fe + r)nf =1 ^>(r,)^' ) (r J -), (A2) 



m^t;}}(r)^^Y. F ^ P ^nt}Kr) = N D f nf =1 dr,j(^r,+r)nf =I |^ 4 (r,)| 2 . (A3) 



Fi has the following properties: it is real and has the dimension of a density; it is normalized to N D , i.e. J drFi(r) = 
N D ; \im r ^± 00 Fi(r) = 0; it is a " multiconvolution" of all single densities, unless the one associated to orbital i. The 
first three points permit to explicit the physical meaning of ■^■^(r): it is the probability that particle i has position 
r, according to the cm. coupling with every other particles and their probability distributions. 

Fp^pi^p is only due to exchange effects. In the following, as explained in £ )II Bl and i jlll B[ we neglect the pure 
exchange effects and thus Fp^pi-tp. We obtain 

f N P 2 

/ dr 1 ...dr N 5(R)ipl lt (r 1 ,...,r N ) V — ^-i/> mt (ri, . . . , r^) (A4) 
J ^ 2m 

h 2 1 N r 



h 2 1 1 



2mNT aux *(Q) R T aM (R) 



R - ^E /*^[{^f}](r)x|^ t (p) 
i=i J 



We now insert this result in E Akin , Eq. (|TU|) . and keep only the real part, i.e. $le(E Akin ), as justified in i]IIIBl It 
only remains the cm. correlations contribution: 

h 2 N f 1 

E Akm E cm = -—Y^Jdr ( —— Fi[{<p$}](T) - l) x ^ t (r)A r ^„t(r) (A5) 

( l l 



2mN \Y aux *{0) T aux (H) 
- Purelmaginary[{tpi nt }] 



where 

ti 2 N f 1 

Purelmaginary[{^ nt }} = i—^mJ2 J dr ___J!i[{^}](r) x ^ t (r)^(r) 

(A6) 

is a pure imaginary functional which counter acts the pure imaginary part of the first line of (|A5[) . 

Ft is interesting to underline the Physics in energy considerations. But it is not a fundamental quantity for the 
potential (that will be obtained by variation of E cm ). We thus introduce a more fundamental quantity, which will 
appear in both the cm. correlations energy and potential, namely the "two particles cm. correlations functional" 
defined in Eq. (J20J) , which is linked to Fi by the relations 

VZ ¥= * : FiiWZuW) = / dr\<p l int (v)f x + O, (A7) 

^l* ^ = ftnt( r ) X /^[{^int IJ( r + r )• 

When (|A7[) is inserted in (| A5|) . we obtain the form (fT5| for the cm. correlation energy, presented in ijlll Bl 
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Appendix B: The multiconvolution theorem. 

We define the Fourier transform T of an integrable function L : 5Re — > ^sm as 



Vr,se5ie: T[L](s) = J dr 



and the inverse Fourier transform T 1 of a function L : SJie — > 5m as 



Vr,s G 3?e : r _1 [^](r) = / ds 



L(r) 



L(s). 



We start from (ivT + 1) integrable functions c/i : Sfte — > 3to and define the "multiconvolution" 

if 

C'Kfl'JKr) = / *i- • -dr K .gi(ri) x ■ ■ • x qk^k) x 9k+i(~ ^2 r t - ?)• 

^ i=l 

We calculate the Fourier transform of C: 

T[C](s) = / dr e- 2 — f C[{ 5i }](r) 



A" 



dr dn . . . dv K e 27ris r flri(ri) x • ■ ■ x qk^k) x ffir+i(- X] r * - f ) 

i=i 

/" dwdv x ...dv K e - 2 »'(-££i r --) 5l ( ri ) x ••• x ffjf ( rif ) x 5A - +1 (w) 
dr 4 e 2 " s r ' ffl (r0 



(Bl) 



(B2) 



where the second line has been obtained doing the variable change w = — X^i r i~ *S which implies dw = — df . 

This is the "convolution theorem" [36j generalized to multiconvolutions, which states that the Fourier transform 
of a multiconvolution is the product of the Fourier transforms of each functions that enter into the multiconvolution. 
This theorem permits to fasten the numerical calculation of C[{<7j}](f). Indeed, it translates a multiconvolution to a 
multiplication of functions. Once every T\gi] are calculated, C[{g,}] is given by the inverse Fourier transform of their 
direct product: 

c[{9i}]{v) = T- 1 [n^T^K-s)! ( f ) = T' 1 fe 1 ™! (-?)• 



Note that this relationship is only valid for the form (|B1[) of the Fourier transform. For forms normalized in other 

1 K 

ways, a constant scaling factor, typically proportional to (2tv) d ^ 2 , will appear. 



Appendix C: Some properties of Y aux when N becomes very large. 

The limit where the cm. correlations become negligible is obtained when TV becomes very large, as mentioned in 
£ ]III CI Indeed, the "two particles cm. correlations functional" /t,/^i then tends to become constant and delocalized 
in the whole space. We note 

lim fi i-a — Constant. (CI) 

The normalization condition ([24| thus implies, when N is very large (remind that Y aux is implicitely dependent of 
N) 

lim |r— (0)| 2 = lim / dr dr' |^„*(r)| 2 |^(r / )| a /i,^i[{^ i,I }](r + O = Constant. (C2) 
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When those results are inserted in E cm , Eq. (TIB)) , we see that its first and obviously third lines become null, and that 
its second line becomes proportional to x ra „ Ar rau l^ x |r au;r (0)| 2 , which must tend to zero so that 



R=0 

E cm can be neglected, as it should be when N becomes very large. This implies 



1 1 A 1 

lim — x - — Ar- 



n^+ooN T aux *(0) T aux (R) 



x \T aux (0)\ 2 -> 0. (C3) 

R=0 



To obtain a second relation, we note 1Z the region of space where the system has a non zero density and V the 
corresponding volume. For very large systems, we have 

bintWI 2 for reft, 

w 0, for r i U. (C4) 

Inserting those results in the definition (j2"3")l of /y^i gives 

fidM r ) ~ for r e ft, 

fa 0, for r <^ ft. (C5) 

In the general case, we have V < fciV, where fc is a constant (it is for saturating systems, as nuclear ones 0, that V 
becomes close, but still inferior, to kN). Thus 

lim /i,;^i(f) = +oo, (C6) 

N— »+oo 

whatever the dimension in which the calculation is done (but limAr_ ! . +00 jji? fi,i^i(r) — > \impf^. +00 y, which tends to 
zero). As a consequence of Eqs. (|C1|) . (IC2[) and (|C6[) . we deduce a second relation that r aux should satisfy: 

lim \T aux (0)\ 2 -> +oo. (C7) 

iV— >-+oo 
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